(* Коэффициенты лин-й модели ЛА *)
a1={0.635,1.38,2.5};
a2={5.47,25.5,54.4};
a3={2.72,9.2,10.6};
a4={0.269,0.527,0.905};
a5={3.26,13.1,10.4};
a6={0.709,1.13,1.6};
a7={0.043,0.0709,0.0554};
b1={3.1,7.32,12.6};
b2={20.2,61.7,176};
b3={17.6,51.2,33.5};
b4={0.0719,0.036,0.024};
b6={0.0571,0.126,0.178};
b7={0.0649,0.0134,0.0016};

(* Выбор режима 1..3 *)
rezh=2;

(* Передаточные числа АП крена с ЖОС *)
tr=3;(* зададим время регулирования *)
myu=(18-b1[[rezh]]*tr)/(b3[[rezh]]*tr)
ii=108/(b3[[rezh]]*tr*tr)
nyu=216/(b3[[rezh]]*tr*tr*tr)

(* Матрмца состояния свободного самолета *)
MatrixForm[MurKrap0={{(s+b1[[rezh]]),0,b3[[rezh]]},{-1,0,0},{myu,nyu/s,-1}}]
(* Матрица с заданным углом крена *)
MatrixForm[MurKrap={{(s+b1[[rezh]]),0,b3[[rezh]]},{-1,s,0},{myu,(ii+nyu/s),-1}}]

(* Расчет передаточной функции *)
w=Simplify[Det[MurKrap0]/Det[MurKrap]]

Wt=Det[MurKrap0]/Det[MurKrap]

(* Перейдем от s к комплексной форме *)
wk=w/.s->I*omega

f=(Abs[wk])^2

Plot[f,{omega,0,5},PlotRange->{0,1}]

p=NIntegrate[f,{omega,0,Infinity}]

NIntegrate[f,{omega,0,∞}]

omega0=0;
shag=0.1;
While[shag≥10^-5,{While[(NIntegrate[f,{omega,0,omega0}]-0.95*p)<0,omega0=omega0+shag],omega0=omega0-shag,shag=shag/10}]
omega0

p95=NIntegrate[f,{omega,0,omega0}]

(* ПЕРИОД КВАНТОВАНИЯ *)
Tkvant=Pi/omega0

ht=InverseLaplaceTransform[Wt/s,s,t]

Plot[{ht,1,0.95,1.05},{t,0,10},PlotRange->{0,1.1}]

hf[t_]=Abs[ht]

(* Время регулирования и перерегулирование *)
k=10000;(*кол-во промежутков, где находим решение*)
dt=10/k;
hust=hf[10000](*установившееся значение*);
hmax=hf[0](*зададим начальное макс. значение*);
(*Численное решение уравнений h(t)=0.95 и h(t)=1.05, объединенных в одно |h(t)-1|-0.05=0 и определение максимума*)
For[i=1,i<k+1,{If[Sign[Abs[hf[dt*(i-1)]-hust]-0.05]≠Sign[Abs[hf[dt*(i)]-hust]-0.05],Treg=dt*(i-0.5)],If[hmax<hf[dt*i],hmax=hf[dt*i]]};i++]
sigmaper=Round[100*(hmax-hust)/hust,0.01]; (*перерегулирование*)
Print["Treg = ",Treg," с"] 
Print["Sigma = ",sigmaper," %"]

(*Кол-во периодов квантования за время регулирования*)
Nkvant=Ceiling[Treg/Tkvant]

Stextrop=3(*Степень полинома экстраполяции*);
Masotklon={}(*Массив отклонений*);
hizm=Table[0,{Stextrop}](*Массив последних измеренных (реальных) значений используемых при экстраполяции*);

For[i=Stextrop-1,i<=Nkvant,i++,
{Ttek=Tkvant*i;
hTtek=hf[Ttek];
hizm[[1]]=hf[Ttek-Tkvant*2];
hizm[[2]]=hf[Ttek-Tkvant*1];
hizm[[3]]=hf[Ttek-Tkvant*0];
(*Экстраполяция Метод Адамса*)
hextrop=hizm[[3]]+Tkvant*(23*hizm[[3]]-16*hizm[[2]]+5*hizm[[1]])/12;
Masotklon=Append[Masotklon,Abs[hf[Ttek+Tkvant]-hextrop]];
}];

(*Математическое ожидание отклонения истинного угла крена от экстраполированного значения*)
Matozh=Mean[Masotklon]

(*Среднее квадратическое отклонение отклонения*)
sigma=StandardDeviation[Masotklon]

(* ПОРОГ СРАБАТЫВАНИЯ *)
b=3 (*1..5 - доверительная граница*);
Porog=Matozh+b*sigma